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Abstract 

In this work we investigate the process of pattern formation in a two 
dimensional domain for a reaction-diffusion system with nonlinear dif- 
fusion terms and the competitive Lotka-Volterra kinetics. The linear 
stability analysis shows that cross-diffusion, through Turing bifurca- 
tion, is the key mechanism for the formation of spatial patterns . We 
show that the bifurcation can be regular, degenerate non-resonant and 
resonant. We use multiple scales expansions to derive the amplitude 
equations appropriate for each case and show that the system sup- 
ports patterns like rolls, squares, mixed-mode patterns, supersquares, 
hexagonal patterns. 

1 Introduction 

The aim of this paper is to study the pattern formation for the reaction- 
diffusion system: 



— =V • Ji + r«(/ii - 7nu - 7i 2 u), 

— =V • J 2 + Tv(fl 2 ~ 721U - 722^)- 

at 

where the fluxes Jj have the following nonlinear expressions: 
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Ji =V(u(ci + a\U + bv)), 
J 2 =V(w(c 2 + a 2 v + b 2 u)). 



(1.1b) 



Here u(x, t) and t> (x, i) are the population densities of two competing 
species, and x G Q with = [0,L X ] x [0,L y ]. The above system is supple- 
mented with initial data and following Neumann boundary conditions: 

n • Ji = n • J 2 = when x £ <9$7 . 

The nonlinear diffusion terms describe the tendency of the species to 
diffuse, when densities are high, faster than predicted by the usual linear 
diffusion towards lower density areas. The parameters and c% are the 
self-diffusion and the linear diffusion coefficients respectively, while the pa- 
rameters b and 62 are the cross-diffusion coefficients. For two competing 
species it is natural to suppose all these parameters to be non negative. 

The constants 7^ > are the competitive interaction coefficients, the 
constants fa are the rates at which each species would grow in absence of 
competition, while the parameter V regulates the size of the spatial domain 
(or the relative strength of reaction terms). 

Since the seminal paper of Turing [45) . reaction-diffusion equations are 
one of the best-known theoretical models explaining self-regulated pattern 
formation in many different areas of physics, chemistry, biology, geology, 
etc. Turing showed that the interplay of diffusion and kinetics can destabi- 
lize the uniform steady state and generate stable, stationary concentration 
patterns. However, if the model has a trivializing kinetics, as it is the case 
for competitive Lotka-Volterra kinetics, classical diffusion is not sufficient 
to destabilize the equilibria, no matter what the diffusion rates are, and 
no pattern formation can be observed. Thus, in order to model segrega- 
tion keeping a simple form for the kinetic term, Shigesada, Kawasaki and 
Teramoto [43J proposed the nonlinear evolution system (II. lj) . 

Strongly coupled parabolic systems with nonlinear diffusion terms of the 
form given in (jl.ip have been used to model different physical phenomena 
and appeared in many context like chemotaxis jHTJ [2], ecology [271 H5 EH 
EH [Ml EH SO] , social systems [351 EH IS] , turbulent transport in plasmas 
|17j , drift-diffusion in semiconductors [111 [TJ [16] , granular materials [H [23] 
and cell division in tumor growth [42] . 

The system (jl.ip has been extensively investigated from the mathemat- 
ical point of view. In |9l [TO], H7j global existence and regularity results have 
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been obtained, while in [32] the existence of non constant steady state solu- 
tions in the time independent case was investigated. The proof of existence 
and stability of traveling wave solutions has been obtained in [49J. The ex- 
istence of positive steady-state solutions in relation to large cross- diffusion 
coefficients has been discussed in [3§] . Some families of exact solutions have 
been constructed in [12] using the Lie symmetry approach while Lyapunov 
functionals have been used in [20, 36J to obtain stability and instability cri- 
teria of the zero solutions of cross-diffusion systems. Prom the numerical 
viewpoint different numerical schemes have been proposed to solve reaction- 
diffusion systems with nonlinear diffusion, see [31 [22j [251 E EE]- 

The importance of the cross-diffusion, relatively to pattern formation, is 
extensively discussed in [H [36] from both the experimental and the theoret- 
ical point of view. In these papers the authors report many experiments of 
interest to chemists where cross-diffusion effects can be significant: they ob- 
tain the minimal conditions for pattern formation in the presence of linear 
cross-diffusion terms, demonstrating that relatively small values of cross- 
diffusion parameters can lead to spatiotemporal pattern formation provided 
that the kinetics is sufficiently nonlinear. 

The focus of this work is to describe the mechanisms of pattern formation 
for the system with homogeneous Neumann boundary conditions in 

a 2D domain. The crucial difference with the ID case analyzed in [26] 
consists in the possibility that bifurcation occurs via a simple or a multiple 
eigenvalue, see below and [33 \ I38| , [T5] . We shall perform a weakly nonlinear 
analysis close to the bifurcation state using the multiple scales analysis: this 
will give the equations which rule the evolution of the pattern amplitude 
near the threshold. A systematic approach to derive the normal forms and 
the correspondent amplitude equations for flows at local bifurcations can be 
found in [T31 [3H [29j [131 HI]- A comparison between different methods of 
weakly nonlinear analysis for several prototype reaction-diffusion equations 
is given in |48j . 

The paper is organized as follows: in Section [2] a linear stability anal- 
ysis of the system (jl.ip is performed and the cross-diffusion is proved to 
be responsible for the initiation of spatial patterns. In Section [3] the am- 
plitude and the form of the pattern close to the bifurcation threshold are 
investigated by using a weakly nonlinear multiple scales analysis. In par- 
ticular, when the homogeneous steady state bifurcates to spatial pattern 
at a simple eigenvalue, we derive the cubic and the quintic Stuart-Landau 
equation which rules the evolution of the amplitude of the most unstable 
mode in the supercritical and subcritical case respectively. In these cases 
the system supports patterns such rolls and squares. On the other hand, 



3 



when the bifurcation occurs via a double eigenvalue more complex patterns 
arise due to the interaction of different modes (for this reason they are called 
mixed mode patterns). The corresponding evolution systems for the ampli- 
tudes of the pattern are obtained and analyzed. A particular type of mixed 
mode patterns are the hexagonal patterns, which arise when a resonance 
condition holds. The evolution system for the amplitudes of the hexagonal 
patterns is proved to show bi-stability and the phenomenon of hysteresis can 
be observed. 

In all the considered cases the solutions predicted by the weakly nonlinear 
analysis are compared with the numerical solutions of the original system. 
Close to the threshold they show a good agreement. 



2 Cross-diffusion driven instability 

In this section we shall investigate the possibility of pattern appearance for 
the system (11. ip . In Subsection 12.11 we shall determine the critical value for 
the bifurcation parameter and the critical wavenumber via linear stability 
analysis. This will be done ignoring the geometry of the domain and the 
role played by the boundary conditions. 

This role will be considered in Subsection 12.21 where we shall obtain the 
range of the unstable wavenumbers of allowable patterns strictly depending 
on the domain geometry. Since the degeneracy phenomenon can occur, the 
situation is much more involved than in the 1-D domain treated in [261 . 



2.1 Main results on the destabilization mechanism 

In order to stress the role played by the cross diffusion term in the pattern 
forming process, the kinetics is chosen of the simplest form, namely the 
competitive Lotka-Volterra, i.e. all 7« > 0. We shall only analyze the 
coexistence equilibrium: 

/ v ( W722 - ^2712 M2711 - M1721 \ / 01 n 

(uo,v Q )=\ , , (2.1) 

V 711722 - 712721 711722 - 712721 / 

as this is the only steady state relevant for pattern formation. Therefore in 
this paper we shall assume the following conditions: 

M1722 - Ai27i2 > 0, ^2711 - M1721 > 0, 7H722 - 712721 > . (2.2) 

The third condition (weak interspecific competition) is necessary for the 
stability of the equilibrium (12. ip . 
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Upon linearization of the system (jl.ip in a neighborhood of (uq,vq), 
namely: 



W = FKw + L>V 2 w, where w = ( U U ° ) , (2.3) 



v-v 



and where: 



if = " 711U0 " 7l2no , (2.4) 

V -721^0 -722^0 J 

D _ _ ( ci + 2ain + 6-u ou \ / 2 ^ 

\ b2Vo C2 + 2a 2 vo + 62^0 /' 

we look for solutions in the form e lk ' x+<T *. Substitution in (|2.3p leads to the 
following dispersion relation, which gives the eigenvalue a as a function of 
the wavenumber k = |k|: 

a 2 - g(k 2 )a + h(k 2 ) = , (2.6) 

where 

g (k 2 ) = k 2 tv(D) - T tr{K) , h(k 2 ) = det{D)k A + Tqk 2 + r 2 det(i-Q , 

(2.7) 

and 

q = 7iiUo(2a 2 i>o + c 2 ) + 722V (2aiii + c%) + 6^0(722^0 - 721^0) , oN 
+ »2^o(7n u o - 712V0) • 

For the Turing instability to be realized and spatial patterns to form, Re(cr) 
must be greater than zero for some k ^ 0. Since the polynomial g(k 2 ) > 
0, \fk 7^ (in fact ti(K) < being (uo,vq) stable and tr(D) > 0), Re(a) 
will be positive for that / at which h(k 2 ) < 0. The steady state is 
marginally stable at some k = k c where: 

mm(h(k 2 )) = 0. (2.9) 

As the minimum of k is obtained when: 

k * = (2 ' 10) 

one has to require that q can become negative. In the above formula the 
matrix D bc is the matrix D defined in (12.511 evaluated at b = b c . 
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Prom the expression (|2.8j) of q where it is apparent that the first two 
terms are non negative, it follows that the only potential destabilizing mech- 
anism is the presence of the cross-diffusion terms. Moreover only one of the 
last two terms in (|2.8j) can be negative, due to the conditions on positiveness 
and stability of (uq,vq). Therefore when 722^0 — 721^0 < (verified in the 
hyperbolic sector Pi on the left of Fig. 12. ip . b has a destabilizing effect and 
62 acts as a stabilizer. Alternatively, when 711 uo — 712^0 < (region P 2 
on the left of Fig. I2.ip . b and 62 exchange their role. In the remainder of 
this paper we shall choose the kinetic parameter set in the first case and the 
cross-diffusion coefficient b as the bifurcation parameter. 

Let us now write q = —ab + f3, where the positive quantities a and (3 are 
defined as: 

a = ^0(721^0 - 722^0) , 

P = r yiiu (2a 2 v + c 2 ) + 722^o(2ain + c\) + &2«o (711^0 ~ 7i2^o) • 

Substituting b = f3 /a + £ in the condition for the marginal stability (|2.9p 
leads to the following equation for £: 

a 2 

■ £ - VQ(2a 2 vo+c 2 )i - [vqP I 'a{2a 2 vo + c 2 ) 



4det(#r - u ^" u '^ L ^^v— ■ (2n) 
+(2ain + c\)(2a 2 vo + b 2 u + c 2 )] = 0. 

Then the critical bifurcation value is: 

b c = /3/a + t, (2-12) 

where £ + is the positive root of equation (|2,lip in such a way that q < 0. 
The results of this section can be summarized in the following theorem. 

Theorem 1 Suppose that (uo,Vo), as given in (|2.ip . is a stable equilibrium 
for the competitive kinetics of the system (jl.ip . i.e. assume conditions (|2.2p 
to be valid. 

Moreover assume that 711M0 — 712^0 < 0. If: 

b > b c , 

where b c is defined in (|2.12p . then (uq,vq) is an unstable equilibrium for the 
reaction- diffusion system (jl.ip . 

An analogous statement would hold when 722^0 — 721^0 < 0. 
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Figure 2.1: Left: Pattern forming regions in the plane (7i2>72l)- Right: 
Growth rate of the k-th mode. 



Proposition 2 Suppose that (uq,vq) is a stable equilibrium for the compet- 
itive kinetics of the system (|l.ip . 

If 722^0 ~~ 72i u o < 0, the Turing parameter space of instability is the 
bounded region P\ on the left of Fia X2l\, Alternatively, if 722^0 — 721^0 < 0, 
the Turing parameter space of instability is the bounded region P2 ■ 



2.2 Instability bands and degeneracy 

When the domain is finite, the condition b > b c is not enough to see a 
pattern emerging. In this case, in fact, k c might not be a mode admissible 
for the domain and the boundary condition. However when b > b c , there 
exists a range (kf,^) of unstable wavenumbers that make h(k 2 ) < and, 
correspondingly, Re(a) > 0, see on the right of Fig J2.1i It is easy to see that 
the extremes of the interval of unstable wavenumbers, k\ and k%, where 
h(k 2 ) = 0, are proportional to T. It follows that F must be big enough to 
find at least one of the modes allowed by the Neumann boundary conditions 
within the interval [kf, fe|]. 

In a rectangular domain defined by < x < L x and < y < L y , the 
solutions to the linear system (12. 3ft with Neumann boundary conditions are: 

fmn e CT ( fc ™)* cos — x cos — y ), (2.13) 
,«=n \ L x J \ L y J 



w = 

m,n£N 



/ m7r\ 2 / mr\ 



where f m>n are the Fourier coefficients of the initial conditions; the values 
a (hmn) are derived from the dispersion relation (12. 6p . The occurrence of a 
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pattern emerging as t increases, therefore depends on the existence of mode 
pairs (m,n) such that: 



i.e. for b > b c and T sufficiently large. In what follows we shall restrict 
ourselves to the case when there is only one unstable eigenvalue, admissible 
for the Neumann boundary conditions, that falls within the band (k\, fe) m 
the sense of Eq. (|2.15p . We shall denote this admissible eigenvalue with k c 
to distinguish from the critical value k c . 

In a two dimensional domain, g 1VGI1 Ivq £ [^1 5 ^2 

], one, two or more pairs 

(m, n) may exist such that the condition 



is satisfied and in this case the eigenvalue a will have single, double or higher 
multiplicity respectively. The multiplicity of the eigenvalue, and therefore 
the type of linear patterns we could expect, strictly depends on the dimen- 
sions of the domain L x and L y . 

In Fig l2,2l we show the pattern which forms starting from an initial datum 
which is a random periodic perturbation about the steady state (uq,vq). 
For the parameters we have picked for this simulation, one has (uq,vo) ~ 
(1.67, 0.92) while the critical value of the bifurcation parameter is b c = 7.192. 
In the rectangular domain with L x = \^2tt and L y = 2ir only the mode 
k% = 6 is admitted by the boundary conditions. The eigenvalue predicted 
by the linear analysis is single: in fact there only exists the pair (2, 4) which 
satisfies the condition (|2.14p . 

All the numerical simulations showed in this paper are performed by 
using spectral methods. Here we have employed 32 modes both in the x 
and in the y axis. However the use of a higher number of modes in the 
scheme (we tested the method up to 128 modes) does not appreciably affect 
the results. Notice that in all the figures representing the spectrum of the 
solutions, for a better presentation of the results, the amplitude of the zero 
mode, corresponding to the equilibrium solution, has been set equal to zero. 

3 Weakly nonlinear analysis 

In this section a weakly nonlinear analysis is carried out to obtain the am- 
plitude equations describing the dynamics near the critical bifurcation state 





(2.17) 
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Figure 2.2: Left: The species u. Right: The species v. The parameters are 
Ml = 1-2, M2 = 1, 7n = 0.5, 7i 2 = 0.4, 721 = 0.38, 722 = 0.4, a x = 0.01, 
a 2 = 0.001, ci = 0.1, c 2 = 0.2, T = 28.05, b 2 = 1.1, 6 = 7.264 > 6 C = 7.192. 




Figure 2.3: Spectrum of the numerical solutions in Fig. 12.21 
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. The method of multiple scales (as introduced in [37]), together with an 
asymptotic analysis of the system close to its marginal stability are employed 
to determine the near-critical bifurcation structure of the patterns [13] • 

Near the threshold the amplitude of the pattern evolves on a slow tempo- 
ral scale, and therefore we shall introduce new scaled coordinates separating 
the fast time t and slow time e.g. T = et; here the control parameter e 
will measure the distance of the system from bifurcation, see (I3.7P below. 
The solution of the original system (jl.ip is written as an expansion in e 
and the leading order term of this expansion is shown to be the product of 
the basic pattern (the critical solution of the linearized system (|2.3p ) and a 
slowly varying amplitude (see [28l [29] ) . We shall confine ourselves to pat- 
terns which are modulated in time but not in space, so that in our analysis 
we do not take into account the slow spatial scale. 

If one defines the linear operator: 

C b = TK + D b V 2 , (3.1) 

where K and D b are given in (|2.4p and (|2.5p . and if one introduces the 
following bilinear operators acting on (x, y) with x = (x u ,x v ) and y = 

(y u ,y v )- 

n , \ r ( -2 7 n* V - 7i2 (x u y v + x v y u ) \ ,o ^ 

U K [*,y) - 1 ^ _2 722 xV-72l(xV + zV) )' 1 ' 

6 / 2 ai x u y u + b(x u y v + x v y u ) \ 

~&D\x-,y) y 2a 2 x v y v + b 2 {x u y v + x v y u ) J ' { J 

the original system (jl.ip can be rewritten separating the linear and the 
nonlinear part as follows: 

a t w = £ b w + ^Qx(w, w) + Iv 2 Q^(w, w), (3.4) 



with w defined as in (|2.3j) . 

Let us introduce the multiple time scales: 

Ti T 2 T 3 T 4 

t = — + ^ + + ^ + ■■■ (3.5) 

e e l e 4 e 4 

and expand accordingly both the solution w and the bifurcation parameter 
b: 

w = ewi + e 2 W2 + e 3 W3 + e 4 W4 + 0(e 5 ), (3.6) 

b = b c + Eb^ +e 2 b& +e 3 b® +E 4 b^ + 0(e 5 ). (3.7) 
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Expansion (|3.T|) can be considered as the definition of the smailness parame- 
ter. In the rest of this paper, we shall always measure the distance from the 
threshold using, as unit, the critical value b c . This means that, when dif- 
ferent from zero, we shall set &W = b c . Substituting (|53]l - ([57f]l into the full 
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system (jl.ip . the following sequence of linear equations for w, is obtained: 
0{e) : (3.8) 



£ fcC wi = 0, 



0(e 2 ) : (3.9) 



£ b w 2 = F = ^ - \(Q K + V 2 Q£ J (wlwO 

0(e 3 ) : (3.10) 
£^w 3 = G = ^ + ^-(Q Jf + V 2 Q^) (w l5 w 2 ) 



0(e 4 ) : (3.11) 

<9wi <9w 2 <9w 3 r-72/-i ¥ 



£ f> w 4 = H = ^ + ^ + ^- (Q x + V 2 Q c b ) (wi, w 3 ) 



U\V2 + U 2 Wl 





\ (Qk + V 2 Q^ fec ) (w 2 , w 2 ) - &« V 2 f 
w 3 + W 2 + 

6 (3) V 2 



f ""0 



- 6 (2) V 2 f ^ 

0{e b ) : (3.12) 
foC <9wi chv 2 9w 3 9w4 

Qk + V 2 Qv bc ) (wi, w 4 ) - (Ck + V 2 Qp foc ) (w 2 , w 3 ) 
0° 0° ) ( 5(1)y2w 4 + ^ (2) V 2 w 3 + 6( 3 )V 2 W 2 + fc^VV 

& (l)y2 ^ ^1^3 + U 2 V 2 + U 3 f 1 ^ _ & ( 2 )y 2 ^ ^1U 2 + U 2 Vl 

- 6 (3)V 2 f ^ ' 
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The solution of the linear problem (j3.8[) satisfying the Neumann bound- 
ary conditions is given by: 

m 

Wi = ^2A i (T 1 ,T 2 )pcos(<p i x)cos(il) i y) , (3.13) 
i=l 

where m is the multiplicity of the eigenvalue, A% are the slowly varying 
amplitudes (still arbitrary at this level), while p, which is defined up to a 
constant, is explicitly given by the following formulas: 

, , Ke r( r* - „ = ( i ) , with M = -g^" , 

(3.14) 

where Df-,Kij are the i, j-entries of the matrices Z? fcC and K. 

We shall restrict our analysis to cases where the multiplicity is m = 1 or 

2. 

3.1 Simple eigenvalue, m = 1 

If the eigenvalue is simple, i.e. m = 1, the solution (|3.13|) at the lowest order 
reduces to: 

wi = A(Ti,T 2 )pcos((f>ix) cos(^iy) . (3.15) 

Substituting this expression into the linear equation (|3.9p . the vector F is 
made orthogonal to the kernel of the adjoint of C bc simply by imposing 
Ti = and &W = O.The solution of (USD can therefore be obtained (see 
(??)) and substituted into the linear problem f|3.10|> at order e 3 . The vector 
G, given by (??), contains secular terms, and therefore it does not auto- 
matically satisfy the Fredholm condition. By imposing the compatibility 
condition the following Stuart-Landau equation for the amplitude A(T 2 ) is 
finally obtained: 

^ = aA-LA\ (3.16) 
where the expression of a and L are in terms of the parameters of the system 

In the pattern-forming region the growth rate coefficient a is always 
positive. Therefore one can distinguish two cases for the qualitative dy- 
namics of the Stuart-Landau equation (|3.16p . The supercritical case, when 
L > 0, and the subcritical case, when L < 0. Since the expression for L 
as a function of all the parameters of the original system is quite involved, 
we numerically determine the curves in the space (7125721) across which L 



13 



721 0.2 



0.05 



0.15 



0.25 



0.35 



0.4 



0.3 



0.1 




0.1 



0.2 



0.3 



0.4 



Figure 3.1: Within the pattern forming region, the zones of subcritical (I and 
III) and supercritical (II, dashed) bifurcation are drawn. The parameters 
are fix = 1.2, fi 2 = 1, Til = 0.5, 722 = 0.41, c\ = c 2 = 0.2, a\ = a 2 = 0.1, 
62 = 0.154. 

changes its sign (all the other parameters being fixed). These curves divide 
the pattern forming domain into three regions, as shown in Fig J3.lt region 

II (dashed) corresponds to the supercritical bifurcation, while regions I and 

III correspond to the subcritical case. 

3.1.1 The supercritical case 

When a and L are both positive, the solution of the Stuart-Landau equation 
evolves towards the stable stationary state = \J a/L. Therefore: 

Theorem 3 Assume that: 

1. e 2 = (b—b c )/b c is small enough so that the uniform steady state (no, ^o) 
in (|2,ip is unstable to modes corresponding only to the eigenvalue k c ; 

2. there exists only one couple of integers (m, n) such that: 



3. the Landau coefficient L is greater than zero. 

Then the emerging solution of the reaction- diffusion system (jl.ip is given 

by: 



where A^ is the stable stationary state of the Stuart-Landau equation ([3. 161) . 
and p is given in (I3.14p . 





(3.17) 
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Figure 3.2: Supercritical case. Comparison between the numerical solution 
(on the left) and the weakly nonlinear first order approximation of the solu- 
tion (on the right). The system parameters are chosen as in FigGTTJ at the 
point 712 = 0.4, 721 = 0.38, marked with an asterisk. Moreover, T = 40.75, 
b c = 5.37, e = 0.1, b = (1 + 0.1 2 ) b c = 5.43. 

In the first numerical test we choose the set of parameters in such a way 
that the only unstable mode allowed by the boundary conditions is k 2 = 8 
and a supercritical bifurcation arises (i.e. a and L are both positive). In 
particular the parameters are chosen as in Fig l3.ll at the point marked with 
an asterisk. 

On a square domain with dimensions L x = L y = tt, there exists only 
the pair (m, n) = (2, 2) such that formula (I2.14p is satisfied. In this case 
the system supports square patterns. The first order approximation of the 
solution predicted by the weakly nonlinear analysis reads: 

w = e Ar^ p cos {2x)cos(2y) + 0(e 2 ), (3.18) 

which shows a good agreement with the numerical solution of the full sys- 
tem (II. ip . see Fig. (l3,2j) . We have verified that the error in predicting the 
amplitude is 0(e 2 ). 

3.1.2 The subcritical case 

For certain values of the parameters (see regions I and III in Fig l3.ip . the 
Landau coefficient L has a negative value. In these cases Eq. ()3.16j) is not 
able to capture the amplitude of the pattern. This is a typical situation 
where the transition occurs via a subcritical bifurcation. To predict the 
amplitude of the pattern, one needs to push the weakly nonlinear expansion 
to a higher order (for a general discussion on the relevance of the higher 
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order amplitude expansions in the study of subcritical bifurcations, see the 
recent [6] and references therein). 

Performing the weakly nonlinear analysis up to 0(e 5 ) one obtains the 
following quintic Stuart-Landau equation for the amplitude A: 

A A 

?—=aA-LA 3 + QA 5 . (3.19) 
a±2 

One can summarize the analysis as: 

WNL analysis results in the non degenerate subcritical case 

Assume that the hypotheses (1) and (2) of Theorem^ hold and that 

(3) the Landau coefficient L is negative; 

(4) the coefficient Q is positive. 

Then the emerging solution of the reaction- diffusion system (jl.ip is given 
by: 

w = spAso cos{4>x) cos(^y) + 0(e), (3.20) 
where A^ is a stable stationary state of the quintic Stuart-Landau equation 

It is important to notice that, given that the coefficient Q = 0(e 2 ) while 
a and L are 0(1), the equilibria A^ = O^' 1 ). This means that the emerg- 
ing pattern is an 0(1) perturbation of the equilibrium, which contradicts 
the basic assumption of the perturbation scheme (|3.6p . In the subcriti- 
cal case one should therefore expect significant quantitative discrepancies 
between the full system results and the predictions of the WNL analysis. 
Nevertheless in our simulation we have always found a fairly good quali- 
tative agreement of the approximation (|3.20p with the solution of the full 
system; most importantly we have seen that the bifurcation diagram con- 
structed using (|3.19p is able to predict very well phenomena like bistability 
and hysteresis cycle shown also by the full system. 

In Fig. 13.31 we show a comparison between the numerical solution of the 
system (jl.ip and the weakly nonlinear approximation for the choice of the 
parameters corresponding to the point denoted with a square in Fig i3.ll 
In Table 3.1 the values of the amplitudes of the most excited modes of the 
numerical and the approximated solutions furnished by the weakly nonlinear 
analysis are compared. 

In Fig. 13.41 we show the bifurcation diagram for specific values of the 
parameters: the origin is locally stable for b < b c and, when b = b c , two 
backward-bending branches of unstable fixed points bifurcate from the ori- 
gin. These unstable branches turn around and become stable at some b = b s 
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Figure 3.3: Subcritical case. Comparison between the numerical solution of 
(jl.ip with e ~ 0.32 (Left) and the third order weakly nonlinear approximated 
solution (Right) based on the quintic Stuart Landau equation ()3.19p . The 
parameters are chosen as in the Fig. 13.11 at the point 712 = 0.3, 721 = 0.4, 
denoted with a square. Moreover, T w 32, b c = 5.63, k 2 c = 8, e = 0.1, 
b 6.25. 

so that in the range b s < b < b c two qualitatively different stable states co- 
exist, namely the origin and the large amplitude branches. The existence 
of different stable states for one single value of the parameter allows for the 
possibility of hysteresis as b is varied. In Fig. I3.5l we show a hysteresis cycle 
corresponding to a periodic variation of the bifurcation parameter. Starting 
with a value of the parameter above b c the solution jumps immediately to 
the stable branch corresponding to a pattern whose amplitude is relatively 
insensitive to the size of the bifurcation parameter. Decreasing b below the 
value b c the solution persists on the upper branch and the pattern does not 
disappear. With a further decrease of b below b s the solution jumps to the 
constant steady state. To have the pattern formation one has to increase 
the parameter b above b° . 

3.2 Double eigenvalue and secular terms at 0(e 3 ) 

If the multiplicity in A3. 13f) is m = 2 and the so called no-resonance condition 
holds, namely: 



(pi + <f>j ^ <j)j 



or 



and 



(3.21) 



<t>% - 4>j + <t>i 



or 
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Table 3.1: Subcritical case: Parameters as in Fig l3.31 A relatively good 
agreement is observed. 



Modes 


Numerical solution 


Approximated solution 


cos(2x) cos(2y) 


0.0894 


0.0909 


cos(4x) 


0.0420 


0.0252 


cos(4y) 


0.0420 


0.0252 


cos(4x) cos(4y) 


0.2513 


0.1418 




Figure 3.4: The bifurcation diagram of the quintic Stuart-Landau equation 
(|3.19p . The parameters are chosen as in Fig. 13.31 
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Figure 3.5: Hysteresis cycle. The values of the parameters are the same as 
in Fig. [331 Left. Numerical solutions. Right. The spectrum of the solutions. 
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with i,j = 1,2 and % 7^ j, then the Fredholm alternative is automatically 
satisfied at 0(e 2 ) by imposing T\ = and 6^ = (as in the supercritical 
case with m = 1). Secular terms appear into equation (|3.10p at 0(e 3 ), whose 
solvability condition leads to the following two coupled Landau equations 
for the amplitudes A\ and A 2 : 

<jAi - L X A\ + nxAiA%, (3.22a) 
aA 2 - L 2 Al + Q 2 Aj A 2 . (3.22b) 
Therefore: 

WNL analysis results in the degenerate non resonant super- 
critical case Assume that: 

1. e 2 = (b—b c )/b c is small enough so that the uniform steady state (uo, vo) 
in (|2.ip is unstable to modes corresponding only to the eigenvalue k c ; 

2. there exists two couples of integers (mj,nj), i = 1,2 such that: 

k c = <j)\ + ipl where <pi = ——, ipi = — ; 

3. 4>i an d ipi satisfy the no-resonance condition (|3.21|) ; 

4- (supercriticality) the system (13.22H admits at least one stable equilib- 
rium. 

Then the emerging asymptotic solution of the reaction- diffusion system (|1.1|) 
at the leading order is approximated by: 

w = ep(Ai 00 cos(0ix) cos(^iy) + ^200 cos(0 2 x) cos(^ 2 y)) + 0(e 2 ), 

where (-Aioo, ^200) is a stable stationary state of the system (|3.22 jl . 

The stationary solutions of the equations (|3.22ap - (|3.22bp are given by 
the trivial equilibrium and the following eight points: 

* s (W?) ; (3 - 23) 

p{ ±,±) = ( I ^2 + ^1) / ^i + n 2 ) \ 



dAi 
~aW 2 
dA 2 
dT 2 
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where the first coordinate is the amplitude A% and the second one is A 2 . It 
is straightforward to prove that the trivial equilibrium is always unstable. 

The results of the linear stability analysis of the equilibrium points are 
summarized in the following table: 





Existence 


Stability 


p? 


Li > 


L x + n 2 < o 


p£ 


L 2 > 


L 2 + fii < 



LiL 2 - S7iSl2 < 0, 

L\ + fl 2 < always unstable 

L 2 + fix < 

P 3 (± ' ±) 

LiL 2 - > 0, 

Li+fi 2 >0 Li^i +L 2 ^2 + 2LiL 2 < 

L 2 + Oi > 



From the above table it can be easily seen that when P^ 2 ^' exist stable, 
the equilibria P- , with j = 1,2 are unstable. Therefore, when pattern 
forms, there are two possible asymptotic behaviors of the solution: a mixed 
mode steady state solution arising in correspondence of the stable equilibria 
P 3 (±,±) and a single mode steady state solution when Pp , with j = 1, 2 are 
stable. We also notice that, for a square domain, there is symmetry of the 
modes (i.e. <pi = ip 2 and (fr 2 = and one always has L\ = L 2 and = 0, 2 
(this is obvious for symmetry reasons, and can also be seen by inspection of 
the formulas (??)-(??)); which implies that the conditions for the existence 
and stability of P 1 and P^ coincide. 

In what follows we shall perform one numerical test concerning single 
mode patterns, and three tests concerning the case of mixed modes patterns. 

In our first test we consider the rectangular domain where L x = y/2ir and 
L y = ir and with the choice of the parameters as in the caption of Fig l3.6l 
the unique discrete unstable mode is k 2 c = 9. The mode pairs satisfying 
the condition (|2.17p are (0,3) and (4,1). Moreover the equilibria P 2 ^ and 
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Figure 3.6: Comparison between the numerical solution (on the left) and 
the weakly nonlinear first order approximation of the solution (on the right). 
The system parameters are chosen as follows: T = 52.453, fi\ = 1.2, [i 2 = 1, 
711 = 0.5, 712 = 0.4, 721 = 0.38, 722 = 0.41, c x = c 2 = 0.2, a x = 0.01, 
b c = 7.316, e = 0.1, b = (1 + 0.1 2 )6 C = 7.389, a 2 = 0.1, b 2 = 0.7. 



are unstable and only the steady states P^ are stable. The predicted 
asymptotic solution therefore is: 



W — iS^liooPCOS (2y/2x) cos(y) + 0(e 



(3.25) 



where Ai^ is the nonzero coordinate of one of the points P^ (which equilib- 
rium is reached depends on initial conditions). Our numerical tests starting 
from a random periodic perturbation of the equilibrium show that the so- 
lution evolves to the rectangular pattern predicted in (|3.25l) . Figure 13.61 
shows the agreement (with e = 0.1) between the numerical solution and the 
solution expected on the basis of the weakly nonlinear analysis. 

In the second numerical test we consider a square domain with dimen- 
sions L x = L y = \^2tt and choose the parameter values in such a way that 
only the most unstable mode fc 2 = 13 falls within the band of unstable 
modes. The set of parameters is described in the caption of Fig l3.81 The 
uniform steady state is then linearly unstable to the two mode pairs (1,5) 
and (5, 1). With this choice of the parameters the two single mode steady 
states in (13.230 are unstable and the mixed mode steady states in (|3.24p are 
stable. Therefore the predicted equilibrium solution to first order is: 



w = e ( ^loopcos 



/ X 



COS 



5y \ ( 5x 

V = 2J +A200PCOS {V = 2 



cos 



(a) 



+0(e 2 ), 
(3.26) 
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Figure 3.7: Comparison between the numerical solution (on the left) and 
the weakly nonlinear first order approximation of the solution (on the right). 
Regions where u > 1.735 are shaded in both plots. The system parameters 
are chosen as follows: T = 49.95, fi\ = 1.2, fi 2 = 1, 7n = 0.5, 712 = 0.4, 
721 = 0.38, 722 = 0.41, ci = c 2 = 0.2, a x = 0.1, a 2 = 0.015, b c = 4.354, 
£ = 0.02, b = (1 + 0.02 2 )6 C = 4.356, b 2 = 0.2. 



where Ai 00 ,A2 00 are the coordinates of one of the equilibrium points in 
(|3.24p . Which of these equilibrium points is reached clearly depends on 
the initial conditions. In Fig l3.7l the comparison between the numerical 
simulation of the original system and the expected solution with e = 0.02 
shows a good agreement and the values of the most excited modes of the 
numerical solution and of the solution computed using the weakly nonlinear 
analysis are quite similar, respectively 0.0278 and 0.0282. 

In the third numerical test we consider again the square domain L x = 
L y = \^2tt and the parameter values are chosen as in the caption of Fig J3.8t 
in such a way that only the most unstable mode = 8.5 falls within the 
band of unstable modes. In this domain the uniform steady state is then 
linearly unstable to the two mode pairs (1,4) and (4, 1). As in the previous 
test, the only stable states are the mixed mode steady states in (|3,24p and 
the predicted equilibrium solution, truncated at the first order, is: 



W = £ ( A loo pCOS 



/ X 

W2 



cos [7=2 ) + A2oop cos [7=2 ) cos KT2 ) ) + ° (e2) ' 

(3.27) 

where ^100,^200 are the coordinates of one of the equilibrium points in 
(|3.24p . In Fig l3.8l we show the comparison between the numerical simula- 
tion of the original system and the expected solution with e = 0.1. The 
two solutions are qualitatively similar, in particular the numerical solution 
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Figure 3.8: Comparison between the numerical solution (on the left) and 
the weakly nonlinear first order approximation of the solution (on the right). 
Regions where u > 1.735 are shaded in both plots. The system parameters 
are chosen as follows: T = 24.517, \i\ = 1.2, \ii = 1, 711 = 0.5, 712 = 0.4, 
721 = 0.38, 722 = 0.41, ci = c 2 = 0.2, a x = a 2 = 0.01, b c = 3.024, e = 0.1, 
b = (1 + 0.1 2 )6 C = 3.054, 6 2 = 0.1. 



Table 3.2: 


Modes 


Numerical solution 


Approximated solution 


cos(x/V2) cos(4?//\/2) 


0.0519 


0.0688 


cos(4x/^/2) cos(y/V2) 


0.0519 


0.0688 


cos(3x/V2)cos(3y/V2) 


0.0560 


0.0092 



evolves to a mixed mode steady state. However one can see a significant 
quantitative discrepancy. A closer look, see Table 13. 2^ reveals that this dis- 
crepancy is due to the presence of the subharmonic (3, 3) which is of the 
same order of magnitude of the modes (1,4) and (4,1), while the weakly 
nonlinear analysis predicts to be 0{e 2 ) (see formula (??))• It is interesting 
that the same discrepancy was found in a numerical test performed in [15] for 
a different type of reaction-diffusion system, where the uniform steady state 
was linearly unstable to the same two mode pairs (1,4) and (4, 1). Noticing 
that the subharmonic (3, 3) corresponds to the discrete eigenvalue k 2 = 9 
(according to the formula (j2.15fl ) . one might conjecture that the closeness 
to k 2 = 8.5 makes the linear damping mechanism unable to overcome the 
quadratic excitation mechanism coming from the interaction of the main 
harmonics (1,4) and (4,1). 
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Figure 3.9: The spectrum of the numerical solution on the left of Fig. I3.8L 

In our final test we consider the case of the supersquare as discussed 
in the paper [18], characterized by having k 2 = 5 on a square domain. We 
consider a domain with dimensions L x = L y = 2ir and choose the parameter 
values in such a way that only the most unstable discrete mode k 2 = 5 falls 
within the band of unstable modes allowed by the boundary conditions. The 
uniform steady state is then linearly unstable to the two mode pairs (1,2) 
and (2, 1). With this choice of the parameters the two single mode steady 
states in (|3.23j) are unstable and the mixed mode steady states in (|3.24p are 
stable. Notice that we shall consider only the positive equilibrium P^ + ' + K 
Therefore the predicted solution truncated at the first order is: 

w = e (A loo p cos (z) cos (2y) + A 2oo P cos (2x) cos (y)) + 0(e 2 ), (3.28) 

where ^4ioo 5 ^2oo are the positive values in (|3.24p . In Fig. (|3.10j) we show the 
comparison between the numerical simulation of the original system and the 
solution predicted by the weakly nonlinear analysis with e = 0.05. The two 
solutions are very close and one can see that the accuracy is, as expected, 
0{e 2 )- 

3.3 Double eigenvalue and secular terms at 0(e 2 ) 

If the multiplicity of the eigenvalue is m = 2 and the following resonance 
condition is satisfied: 

<t>i + 4>j = 4>j an d 4>i - ipj = 4>j 

or (3.29) 
<t>i - 4>j = 4>j and ipi + ipj = ipj 
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Figure 3.10: Comparison between the numerical solution (on the left) and 
the weakly nonlinear first order approximation of the solution (on the right). 
Regions where u > 1.735 are shaded in both plots. The system parameters 
are chosen as follows: T = 30.6, fix = 1.2, \i 2 = 1, 711 = 0.5, 712 = 0.4, 
721 = 0.38, 722 = 0.41, ci = c 2 = 0.2, a x = a 2 = 0.1, b 2 = 0.5, b c = 7.181, 
£ = 0.05, b = 7.199. 

with i,j = 1,2 and i 7^ j, then secular terms appear at 0{e 2 ) in (|3.9p . In 
what follows, without loss of generality, we shall thus perform the weakly 
nonlinear analysis imposing that the second condition in (13.29H hold, with 
% = 2 and j = 1. Taking into account these conditions and the relation 
pTfjl . it follows that <p 2 = 2</>i, ip 2 = 0, Vi = \/%<t>i and fa = k c /2, see (??). 
Hexagonal patterns fall within this class of solutions, as shown in [15]. We 
also notice that the above relations imply that L y = \[?>L X . 
The solution (|3. 13[) at the first order reads: 

wi = Ax(Ti, T 2 )p cos(fax) cos(\/30iy) + A 2 (Ti,T 2 )pcos(2fax). (3.30) 

The solvability condition at 0(e 2 ) gives the following system of equations 
for the two amplitudes Ai,A 2 : 

a A x - LAxA 2 , (3.31a) 

oA 2 ~\a\. (3.31b) 

The stationary solutions of the equations (]3.31j) are the trivial equilib- 
rium and Q ± = (± 2a/L,a/L). It is easy to see that the nontrivial sta- 
tionary solutions associated with are always unstable (the eigenvalues of 
the jacobian matrix evaluated at are Ai = — a, X 2 = 2a); therefore the 



dAx 
~dl\ 
dA 2 

dTx 
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weakly nonlinear analysis, at this order, is not able to predict the amplitude 
of the pattern. This is a subcritical bifurcation case and the asymptotic 
analysis has to be pushed to higher order in the amplitude to obtain quali- 
tatively reliable results [6]. To 0(e 3 ) one finds the following system for the 
amplitudes A\ and A 2 : 

dAl - " r a a ,= a3,qaa2 



dT 
dA 2 

dT 



aiAi-LiAi^ + cMi + ZMi^ (3.32a) 
a 2 A 2 -L 2 Al + a 2 Al + p 2 AjA 2 . (3.32b) 



WNL analysis results in the degenerate resonant case Assume 
that: 

1. e = (b — b c )/b c is small enough so that the uniform steady state (uo,vo) 
in (|2.ip is unstable to modes corresponding only to the eigenvalue k c ; 

2. there exists two couples of integers (mj,nj), i = 1, 2 such that: 

72 _ ,2 , ,2 u A - mi7r I - Ui7T 

k c = <t ) i + Wi where <Pi = -j—i yi = ~t~ ' 

3. 4>i and ipi satisfy the resonance condition (|3.29p ." 

4- the system (13.321) admits at least one stable equilibrium. 

Then the emerging asymptotic solution of the reaction- diffusion system (jl.ip 
at the leading order is approximated by: 

w = ep{Ai OQ cos(0ix) cos(V'iy) + A 2oQ cos((j) 2 x) cos(^ 2 y)) + 0(e 2 ), 

where (^4ioo, ^boo) is a stable stationary state of the system (|3.32p . These 
solutions are rolls (when A^ = 0) or hexagons. 

Other than the trivial one, the equilibria of the equations (|3.32|) are the 
points R^ 1 = (0, ±^ —a 2 /a 2 ) which (when exist real) correspond to rolls, 
and the six roots Hf = (Af i ,A 2 i), i = 1,2, 3, of the following system: 

( Al (ata 2 - M2) +A 2 2 (Z^ 2 + L2P1) + 

A 2 (a\a 2 - L\L 2 - j3 2 ai) + £2^1 = 0, (333) 



A\ = ±- {-hAl + UA 2 - a x ) . 
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The roots H i , when exist real, correspond to hexagons. The equilibria 
exist stable when 612 < and L\ > \J —a^jcF^ipx — Pi<J 2 /d<-2)- When these 
points are stable the system (jl.ip support rolls of the form: 

w = e J 4 2 oopcos(2</>ix) + 0(e 2 ), (3.34) 

where A200 is the nonzero coordinate of R^. 

On the other hand a complete analysis of the existence and stability 
of the stationary points H- is too involved to be carried out in general; 
instead we present a detailed numerical study of a typical case. In Fig J3. 11 1 
the parameter set has been chosen such that the only unstable discrete 
mode is k 2 = 4, which corresponds, in the domain L x = 2ir and L y = 
2"v/37r, to the two mode pairs (2, 6) and (4, 0) satisfying the condition (12.14p . 
Numerical simulations, performed choosing as initial condition a random 
periodic perturbation of the equilibrium, shows that the solution evolves 
to the hexagonal pattern on the left of Fig |3.Hl The form of this pattern 
is captured by the following hexagonal pattern, predicted to be a stable 
solution via the weakly nonlinear analysis: 



w 



ep cos(x) cos (V3yj + A 2 \ cos(2x)) + 0(e 2 ), (3.35) 



where Hp = (A^ 1 ,A2i) is a stable equilibrium of the system (|3.32[) . The 
comparison between the numerical solution and the expected solutions f|3.35|) 
is shown in Fig ]3.11[ The WNL analysis predicts very precisely the dominant 
modes, but the overall accuracy is spoiled by the presence of sub-harmonics 
that, as usual for subcritical cases, the WNL analysis underestimates. In 
Fig J3.12l we show the bifurcation diagram. One can see that system (|3.32p 

admits four real roots of the form Hp = ^Ay,A2jj , j = 1,2, where Hp are 

stable and H 2 are unstable. Moreover the system admits as a stable state, 
also the rolls R~ , 

In Fig l3.13l we show the basins of attraction of the three stable points 
R~ , H + and H~ of the system (|3.32p . We have run several numerical tests: 
starting from an initial condition of the form: 

w? = p (a\ cos(^ix) cos(\/3(/>iy) + A° cos(2^i x)) (3.36) 



where (A®, A®) is in the basin of attraction respectively of R~ ,H + or H~ , 
the solution of the original system (jl.ip evolves to the appropriate stable 
solution predicted via weakly nonlinear analysis in each case. The coexis- 
tence of more than one stable points for a single value of the bifurcation 
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Figure 3.11: The numerical solution (on the left) shows hexagonal patterns. 
The expected solution via weakly nonlinear analysis up to the O (e 2 ) is shown 
on the right. The system parameter are chosen: T = 13.6, [i\ = 1.2, fi2 = 1, 
7ii = 0.5, 712 = 0.4, 72i = 0.38 , 722 = 0.39, ci = c 2 = 0.2, a x = 0.01, 
a 2 = 0.001, b 2 = 0.031, b c = 3.599, e = 0.006 and b = 3.62. 




Figure 3.12: The bifurcation diagram of the Stuart-Landau system f)3.32[) . 
The parameters are chosen as in Fig J3.11l 
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Figure 3.13: Basins of attraction of rolls and hexagons. The parameters are 
chosen as in Fig l3.111 

parameter allows the possibility of hysteresis. In Fig l3.14l it is shown this 
phenomenon: starting with a value of the parameter above b c , the solu- 
tion jumps to the stable branch H + and the hexagonal pattern forms (see 
Fig J3.14l with b = 3.67 > b c = 3.599). Decreasing b slightly below the value 
b c , the solution persists on the stable branch and the pattern does not disap- 
pear (see Fig J3.14l with b s < b = 3.52 < b c ). With a further decrease below 
b s = 3.422, the solution jumps to the uniform steady state (see Fig l3.14l with 
b = 3.42 < 6 s ) and persists in this state even after increasing the parameter 
b (see Fig l3.14l with b s < b = 3.51 < b c ). To have the pattern formation 
the parameter must be increased above b c (see Fig J3.14l with b = 3.61 > b c ). 
On the right of Fig l3.14|) we show the spectrum of the solution: the domi- 
nant modes are those prescribed by the weakly nonlinear analysis (i.e. (2, 6) 
and (4,0)). The other modes appearing in the Fourier spectrum are the 
sub-harmonics which are however underestimated by the WNL analysis. 

3.4 Cross-roll instability 

We have also observed the phenomenon of the cross-roll instability. For 
the same parameter set chosen in Fig J3. 11 1 we have performed a numerical 
simulations almost next to the threshold, for e = 10 -4 . Choosing the initial 
conditions as a random periodic perturbation of the equilibrium (where the 
amplitudes of the random modes are small), the solution evolves, as shown in 
Fig J3.15r b). to the roll pattern predicted via the weakly nonlinear analysis: 

w = ey / -a 2 /a 2 pcos{2x) + 0{e 2 ). (3.37) 
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Figure 3.14: Hysteresis cycle. Left. Numerical solutions. Right. The spec- 
trum of the solutions. 
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Roughly at i = 200000 the roll pattern loses its stability and the solution 
evolves to the pattern shown in Fig J3.15r d), 

A new set of rolls, corresponding to the mode (0,7), grows perpendic- 
ularly to the original roll pattern. In the considered domain L x = 2n and 
L y = 2V37T, the mode pair (0, 7) corresponds, in the sense of (|2.14p . to the 
eigenvalue k^ = 49/12 ~ 4.083, which is nearly equal to the predicted most 
unstable eigenvalue k 2 = 4. We shall interpret the formation of the cross- 
roll pattern as due to an intermode competition, following the approach 
presented in [41]. We perform the weakly nonlinear analysis and write the 
solution of the linear problem at 0(e) as: 

wi = p^A(T 1 ,T 2 )cos(2x) + B(T 1 ,T 2 )co S C^vJ) ■ ( 3 - 38 ) 

At 0(e 3 ) the following ODE model that illustrates the nonlinear behavior of 
the amplitudes A and B of the two competing modes k c and k m is obtained: 



A A 

— = o\A — LiA 3 + QiAB 2 , (3.39a) 
d±2 

dB 

— = o 2 B - L 2 B 3 + n 2 A 2 B, (3.39b) 
d±2 

Here we skip all the details as the analysis of this system is similar to the 
case given in Section f3. 2 i For the considered set of parameters, the equilibria 
(±A, 0) and (0, ±i?), each one corresponding to a roll pattern, are unstable. 
The unique stable states are (±A,±B), as shown into the Fig l3.16l which 
corresponds to the cross-roll pattern. 

The above method, however, does not explain the fact that the amplitude 
of the resulting pattern is 0(1) (see Fig l3.15l fd)). which would probably 
require to consider also the competition with the hexagons modes and the 
spatial modulation of the pattern. 



4 Conclusions 

In this paper we have investigated the Turing mechanism induced by non- 
linear cross-diffusion for two coupled reaction-diffusion equations on a two- 
dimensional spatial domain. 

The possibility that the Turing bifurcation occurs via a degenerate eigen- 
value makes the study mathematically involved, but gives rise to a rich 
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Figure 3.15: Cross- Roll instability. Starting from a random initial datum 
the solution evolves to the expected roll pattern (|3.37|) . However a cross-roll 
instability leads to a mixed mode pattern. 
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Figure 3.16: The phase plane of system (|3.39p . The competitive modes are 
fi^? = 4 and = . 

variety of patterns which tessellate the plane and appear as steady state so- 
lutions of the reaction diffusion system. These are rolls, squares and mixed- 
mode patterns, among which there are the supersquares and the hexagons. 

We have obtained the amplitude equations providing a mathematical 
description of the reaction-diffusion system close to the onset of instabil- 
ity. The analysis of the amplitude equations has shown the occurrence of a 
number of different phenomena, including stable subcritical Turing patterns 
or multiple branches of stable solutions leading to hysteresis. In particular, 
hexagonal patterns appear via a subcritical bifurcation and, from its same 
primary bifurcation point, rolls bifurcate supercritically. Therefore, there 
is a region of bistability where both rolls and hexagons are stable; however 
here rolls appear as a transient state due to a spatially modulated cross-roll 
instability that drives the solution toward a mixed modes pattern. Our at- 
tempt at explaining this instability as due to mode competition has been 
partially successful. 

Future analysis might move from stationary Turing patterns to traveling 
patterning waves: when the domain size is large, the pattern is formed 
sequentially and traveling wavefronts are the precursors to patterning. In 
this case the equations governing the amplitude of the spatially modulated 
pattern (one has to consider the slow modulation in space of the pattern 
amplitude) will be the Ginzburg-Landau equation, or systems of Ginzburg- 
Landau equations in the degenerate case; we believe that these systems 
would be of independent and significant mathematical interest. 
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